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A TAYLOR-GALERKIN FINITE ELEMENT ALGORITHM FOR TRANSIENT 
NONLINEAR THERMAL-STRUCTURAL ANALYSIS 

By 

Earl A. Thornton 1 and Pramote Dechaumphai 2 
INTRODUCTION 

The determination of the structural response induced by thermal effects 
is an important factor in many aerospace structural designs. Extreme aero- 
dynamic heating on advanced aerospace vehicles may produce severe thermal 
stresses that can reduce operational performance or even damage structures. 
The performance of laser devices can be degraded by thermal distortions of 
mirror surfaces. The thermal environment in space may cause orbiting struc- 
tures to distort beyond operational tolerances. To predict the structural 
response accurately, effective numerical techniques capable of both thermal 
and structural analyses are required. One technique, the finite element 
method, has been found to be particularly suited for such analyses due to 
its capability to model complex geometry and to perform both thermal and 
structural analyses. 

In the most common approach to determining structural responses induced 
by thermal effects, the thermal-structural analyses are assumed uncoupled, 
and the structural analysis is assumed quasistatic. The uncoupled assump- 
tion means that mechanical deformation terms in the heat transfer energy 
equation are neglected. The quasistatic assumption means that inertia terms 
in the structural equations of motion are neglected. The practical effect 
of these assumptions is that the heat transfer analysis can be performed 
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first, and the resulting temperatures can be used as input to a subsequent 
stress analysis. This approach works well when temperatures change slowly 
as occurs, for example, in an orbiting space structure subject to solar 
heating. Under these circumstances, the uncoupled, quasistatic idealization 
provides an effective approach for finite element thermal-structural analy- 
sis. The principal difficulty encountered in uncoupled analyses is achiev- 
ing compatabi 1 ity between different thermal and structural models. A hier- 
archical finite element approach for integrating uncoupled thermal and 
structural analyses is described in Ref. 1. 

When changes in temperature occur rapidly, the assumptions that justify 
the uncoupled, quasistatic idealization are no longer justifiable. Temper- 
ature changes can occur rapidly due to propagation of thermoelastic waves, 
during vibrations induced by periodic variations of temperature fields, due 
to thermal shocks and in similar circumstances. These types of problems may 
involve resolving wave-like details of the time dependent response for com- 
plex structures. Moreover, if the mechanical coupling terms are retained, 
the equations are inherently nonlinear even in the material's elastic range. 
The finite element method remains a logical solution approach because of its 
capability to represent complex geometries. Tnere is a need, however, for 
effective finite element solution algorithms that can solve large nonlinear 
transient problems efficiently. 

The purpose of this paper is to present a Taylor-Galerkin finite ele- 
ment method for solving large, nonlinear transient thermal-structural prob- 
lems. The method is a thermal-structural application of a Taylor-Galerkin 
algorithm recently developed to solve the conservation equations of invis- 
cid, compressible flow (Ref. 2). In the flow problem, the algorithm is used 
to solve the highly nonlinear Euler equations that includes capturing shock 
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discontinuities in the flow field. Finite element models of flow problems 
usually are quite large with the number of equations typically in the range 
from 3,000 to 30,000 or more. Thus, the algorithm appears to have desirable 
attributes that will make it effective for large nonlinear, transient 
thermal-structural problems. 

The formulation of a nonlinear thermal-structural problem will be pre- 
sented first, then the Taylor-Galerkin algorithm will be described. Next, 
an explicit evaluation of the finite element integrals is described. Final- 
ly, a programming strategy for a vector computer implementation of the 
algorithm is described. 

This research project is supported by NASA/Langley Grant NSG-1321, and 
monitored by Mr. Allan R. Weiting, LAD-Aerothermal Loads Branch. 

THERMAL-STRUCTURAL FORMULATION 

The nonlinear coupled thermal-structural equations for a two-dimension- 
al continuum (ref. 3) are written in the form 

ilJi + . ill . {H} (1) 

at ax ay 

where { u} is the vector of unknowns, { E} and {F} are vectors of stresses 
and heat fluxes, and {h} is a "load" vector. These vectors are given by 
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r 



where u and v are the displacements, p is the density, u and v are 

the velocities, c y is the specific heat and T is the temperature; a^, 

<V t are the stress components; q , q are the heat fluxes; f , f 
y a jr x y 

are body force components per unit volume; &12, 3 22 are coefficients 

that depend on the coefficients of thermal expansion, and Q is the inter- 
nal heat generation rate per unit volume. The first two equations in Eq. 
( 2 ) define the velocity components, the third and fourth equations are tne 
equations of motion, and the fifth equation represents conservation of 

energy. The term in the square brackets in the last line of { H} repre- 

sents the thermal-structural coupling. 

The structural and thermal equations are written in the form of Eq. ( 1 ) 
to resemble the conservation equations of fluid flow. In the fluid context, 
the components of { U} are called the conservation variables, and the com- 
ponents of { E} and { F} are fluxes of mass, momentum and energy across 

the faces of a control volume. In the thermal-structural context, the com- 

ponents of { U} may also be regarded as conservation variables. The stress 
components of { E} and { F} now represent tractions on surfaces of the 

control volume; however, q and q still represent heat fluxes across 

x y 

surfaces of the control volume. 

In formulating the constitutive equations, highly nonlinear relations 
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between stresses, strains and temperatures are permitted as well as nonlin- 
ear relations between heat fluxes, temperatures and temperature gradients. 
Anisotropic materials can be accommodated as well. For simplicity, simple 
constitutive relations for linear, elastic orthotropic materials will be 
presented herein. The stress-strain relations for an orthotropic material 
are expressed as 






_ 


r 




a 

X 


Cn 

Cl 2 

0 



3 U 

3 x 


1 1 

a 

y 

= 

C 2 1 

C 2 2 

0 

< 


3 v 

3 y | 

>' < T -V' 

$22 

T 

xy 


0 

0 

C33 


1 

2- 

^3 U 3 V j 

3 y 3 x 

J 


$12 


where C.. are elastic constants, and T is the reference temperature for 
ij ’o 

zero stress. The heat fluxes are expressed by Fourier's law, 


i 




(3b) 


where k , and k are the thermal conductivities, 
x y 

Equation (1) is solved subject to appropriate initial and boundary 
conditions. The initial conditions consist of specifying the distributions 
of the conservation variables {u} at time zero. The structural boundary 
conditions consist of specifying the displacements or surface tractions at 
all points on the boundary. The thermal boundary conditions consist of 
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specifying temperatures or heat fluxes at all points on the boundary. 
Convective and radiation boundary conditions are incorporated tnrough heat 
fluxes. 


TAYLOR-GALERKIN ALGORITHM 

The solution domain D is divided into an arbitary number of elements 
of r nodes each. Fig. 1 shows typical quadrilateral elements (r=4) used 
in this paper. For simplicity, the finite element formulation will be given 
for a single scalar equation. 


3u + 3E_ + 9F_ = 

at a x ay 


(4) 


where the variables u, E, F and H are analogous to the corresponding vec- 
tor quantities in Eq. (1). Let {u} n denote the element nodal values of 
the variable u(x,y,t) at time t n . The time step At spans two typical 
times t n and t^ + ^ in the transient response. The computation proceeds 

through two time levels, ^x\+l/2 an ^ t n+l* ^ time ^ eve ^ ^+1/2’ va ^ ues 
for u that are constant within each element are computed explicitly. At 
time level t + ^, the constant element values computed at the first time 
level are used to compute nodal values for u. In the time level t n+ ^ com- 
putations, element contributions are assembled to yield the global equations 
for nodal unknowns. The resulting equations are approximately diagonalized 
to yield an explicit algorithm. 

Time Level t n+1 / 2 

To advance the solution to time level t n +i/2» a truncated Taylor 
series yields 
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u(x,y,t n+1/2 ) ■ u(x,y,t n ) ♦ ±- (x,y,t n ) 


(5) 


Then Eq. (4) is introduced on the right hand side of Eq. (5) so that 


u(x ' y>t n+l/2 ) = u(x ' y * t n> 




— (x.y.t ) + — (x.y.t ) 
8 X 3y 


+ I At H(x,y,t) 
2 


( 6 ) 


At time level ^ n +i/ 2 » the dependent variable u ^ x »-y» t n +i/2^ 1S assumec * 
to have a constant value u Q n+1/ ^ 2 within an element. 

At time level t in the response u, E, F and H vary within an ele- 
ment and are interpolated from nodal values. Thus, the following spatial 
approximations are used within an element. 


u ( x »y» t n+ i/ 2 ) = U D 

(7a) 

u(x,y,t n ) = [N(x,y)]{u} n 

(7b) 

E(x,y.t n ) = [N(x,y)]{E} n 

(7c) 

F(x,y,t n ) = [N(x,y)]{F} n 

(7d) 

H(x,y,t n ) = [N(x,y)]{H} n 

(7e) 


where [N(x,y)] denotes element interpolation functions and {u} n is a 
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vector of the element nodal quantities. For a typical quadrilateral element 
the interpolation functions are bilinear, and in a local natural coordinate 
system (see Ref. 4) have the form 

N. =i (1 + £..5) (1 n) i = 1, 2, 3, 4 (8) 

1 4 1 1 

where £.. and n . denote the nodal coordinates (£.j,n^ = 1 1) of the 
element in the £-n plane. The equations for Up n+ ^^ for each element 
are derived by the method of weighted residuals. The spatial approximations 
given in Eq. (7) are introduced into Eq. (6) to give a residual; the result 
is multiplied by a weighting function which in this case is unity. Finally, 
the weighted residual is integrated over the area A of the element. The 
result is. 


n+1/2 

K, - J A [NJ dA {u} n (9) 

- — / A [— ] dA { Ep - c— ] dA {Fp 

2 3x 2 a 3y 

+ jt/ A [N] dA {Hp . 

With Eq . (9), the dependent variable uP + ^ for each element can be com- 
puted explicitly using nodal values for u, E, F and H from the previous time 
t n> The constant element variable u Q n+ ^ may be interpreted as a 
weighted average of an element's nodal values at time t p . A later section 
will discuss the evaluation of the three integrals that appear on the right 
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hand side of Eq. (9). 

In advancing the solution to the next time level, the values of the 
dependent variables on the surface may be required also. Let u $ 
denote the surface value on a typical element edge IJ on the surface, S 3 
in Fig. 1. Following the approach used previously, u s n+ ^ is assumed 
constant on edge IJ at time t n+ i/ 2 ’ time t n , u, E and F vary 

along the edge. Thus the following approximations are used on an element 
edge on S 3 , 


, , . n+1/2 

u(x,y,t n+ i /2 ) = u s 

(10a) 

u(*.y.t n ) " [N(s) j{u) n 

(10b) 

E(x,y.t„) - [N(s)]{E} n 

(10c) 

F(x,y.t n ) - [N(s)][F} n 

(10d) 

H(x,y,t n ) > [t»(s)](H} n 

(10e) 


where [N(s)] denotes the interpolation functions along an element edge. 
Using the method of weighted residuals, the values for u s n ^ are derived 
by integrating over the length L of an element edge. Hence 

L u n+1/2 = /, [N] ds{ u} n - £!/, [N] ds {E} n 
s L 2 

- [LH] ds {F} n + — /, [N] ds { H } n . (11) 

2 L 3y 2 L 
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Thus, Eqs. (9) and (11) can be used to advance explicitly the element and 
surface values of the dependent variables to t n+ i/ 2 * Beginning with nodal 
values of {u} n , { E} n , {F} n and {H} n , Eq. (9) is used to compute constant 
values Up l1+ 1 / ^ for each element. In a similar way, Eq. (11) is used to 
compute constant surface values u s n+ ^ for element edges on boundaries 
that require these values at t n+ y 2 . For example, u $ n+ ^ are required 
on element edges with specified stresses or heat fluxes. These values are 
computed explicitly by looping through all elements and appropriate element 
edges. 

Time Level t p+ ^ 

To advance the solution to t , forward and backward truncated 

n+1 

Taylor series expansions at t n+ jy 2 are usec * t0 wr i te The approximation 

u(x,y,t n+1 ) = u (x,y,t n ) + At — (x,y,t n+1 / 2 ) . (12) 

3t 

Then, following the approach used previously, Eq. (4) is introduced on the 
right hand side to yield 


u(x,y,t n+1 ) = u(x,y,t n ) 

" At ^ x ’ y,t n+l/2^ + 7 ~ ^ x * y,t n+l/2^ * 

dx o y 

+ At H(x,y,t n+1/2 ) . 


(13) 
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If the flux gradients 3E/3x and 3F/3y are evaluated using values from 
the midstep, the gradients will be zero because E and F are constant 
within the element at midstep. To avoid this undesirable situation, the 
flux gradients are interpolated from their values at t n and t n+ ^. Thus 

“ ^ x,y,t n+l/2^ = 6 ~ ( x ’ y ’ t n ) + (1 ' 0) ~ (x ’ y ’ t n+l ) 

~ ( x » y »Vl/2 } = 6 “ < x ’ y >V + (1_0) ~ (14) 

where the interpolation parameter 9 varies from zero to one. The equa- 
tions for the nodal values of { u} can next be derived by the method 
of weighted residuals in the standard way (ref. 4) using the interpolation 
functions as weighting functions. In the process, Eq. (14) is substi- 
tuted into Eq. (13), and the terms containing derivatives are integrated by 
parts. These operations yield the equations for the nodal values of a 
single element, 

[M] {u} n+1 = [M] { u} n + 6 { R x } n +6 { R 2 } n 

+ (1-0) { R 3 } n+1 + (1-0) {R 4 } n+1 + {R 5 }" +1/2 US) 
where [M] = { N} [N] dA (16a) 
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(16b) 


{Rl!" ■ n f. {— f [N] dA (E( n * At J, (!!!) [N] dA (F) n 
a 3X 3y 


<R 2 )" 

= - At / 

s { N) [N] ds (t {£)" 

+ m (F} n $ ) 

(16c) 

( R 3 } n+1 

= At J A 

(— } [N] dA (E) n+1 

+ At /. ft} [N] dA { F} n+1 

( 16d) 


3X 

3y 


(M n+1 

= - At / 

s {N} [N] ds U {E} 

n+1 , r r1 n+1 > 

s + ™(F} S ) 

(l6e) 

(Rs)" +1/2 

= At / A 

{N} dA H n+1/2 

• 

(16f) 


In Eqs. (16c) and (16e) the coefficients a and m are the components of a 
unit vector normal to the boundary. Following usual finite element proce- 
dures, the element matrices given in Eq. (16) can be assembled to form sys- 
tem (global) equations. 

The matrix [M] defined by Eq. (16a) is the element consistent mass 
matrix. The term, consistent, is used to distinguish these matrices from 
diagonal mass matrices that arise from other discretization methods. A 
consistent mass matrix has off-diagonal terms that couple the element vari- 
ables on the left hand side of Eq. (15). The algorithm with consistent mass 
matrices and arbitrary 6 is an implicit scheme. If the mass matrices are 
diagonalized and 0 is taken as 1, the algorithm becomes an explicit 

scheme. 

EXPLICIT EVALUATION OF ELEMENT INTEGRALS 
Element integrals for the Taylor-Galerkin algorithm shown in Eqs. (9), 
(11), and (16) were evaluated in closed form in Ref. 2 to avoid expensive 
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numerical integrations that are customary for quadrilateral elements. For 
two dimensional problems, these element integrals are either in the form of 
integration over the element area or along an element edge. As an example, 
the element mass matrix for a quadrilateral element, Eq. (15), is given by 

[M] = p /l {N} [N] |J | d5 dn (17) 

where { N} is the vector of the element interpolation functions in terms of 
the element natural coordinates K and n. The determinant of the Jacobian 
| J | that appears in the above equation represents the transformation from 
the element global coordinates x-y to natural coordinates S-n. The 

transformation permits the element integration to be evaluated over a 
square. 

The approach developed in Ref. 2 was also used for the evaluation of 
three dimensional hexahedral element integrals. Typical element integrals 
derived using the Taylor-Galerkin algorithm are in the same form as obtained 
for the quadrilateral element where the element integrations are either 
performed over the elanent volume or the element surface areas. As an exam- 
ple, the mass matrix for a hexahedral element is given by 

CM] - J y {N} [N] dV 

= /} P /J { N} [H] |J| <K dn d? (18) 

where C, n, c are the element natural coordinates. 

The CPU time used by the closed form solution of the mass matrix has 
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been investigated and compared with CPU times required for different orders 
of Gauss numerical integration for quadrilateral and hexahedral elements. 
Although CPU time savings for quadrilateral elements are large (a factor of 
50 for all integrals), the CPU time savings for hexahedral elements is even 
more significant. Typical CPU times for evaluation of a hexahedral element 
mass matrix are compared in Fig. 2. The figure shows that the closed form 
solution reduces CPU time significantly for an element mass matrix. Time 
savings in excess of an order of magnitude are obtained from the closed form 
solution in comparison with the popular 8 points Gauss integration method. 
Computational savings from explicit evaluation of other hexahedral element 
matrices are shown in Fig. 3. 

The time savings gained by explicit evaluation of the finite element 
integrals is an important step toward developing efficient finite element 
computations for large three-dimensional problems. 

VECTOR PROGRAMMING STRATEGIES 

The Taylor-Galerkin algorithm used in Ref. 2 was implemented with vec- 
torization strategies specifically for the Langley VPS 32 (a Cyber 205 with 
16 million words of central memory). This computer achieves high computa- 
tional speed when performing operations on long vectors. Vector lengths of 
at least 60 are required to justify vectorization efforts with maximum pay- 
off achieved for vector lengths of 1000 or more. The predominant vector 
lengths in the vectorization scheme are the number of elements in the finite 
element model with occasional operations using vector lengths equal to the 
nunber of nodes. 

The critical vectorization tasks are for those operations that are 
repetitive and performed at every time step. For finite element algorithms, 
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Comparison of execution times for evaluation of 
hexahedral element mass matrix. 



Fig. 3 Comparison of execution times for evaluation of 
hexahedral element matrices. 


these operations are: (1) assembly of element contributions into the global 
system of equations, (2) solution of the global system of equations, and (3) 
application of boundary conditions. 

The assembly of element contributions into the global system of equa- 
tions is the process in finite elements which differs most from other num- 
erical techniques and requires special routines for vectorization. Nodal 
unknowns are stored in one dimensional arrays from 1 to the number of nodes 
in the model, and in general, node nunbering may be arbitary throughout the 
mesh. Assembly of element contributions is performed using the VPS 32 
FORTRAN-suppl ied scatter routine which places an element contribution into 
the proper location in the system of equations based on the element connect- 
ivity. Every element that contains a particular node in its connectivity 
provides its own contribution to the system equations, therefore, the assem- 
bly is an additive operation. "Scattering" alone would merely overwrite a 
previous element contribution. The special vector routine, then does an 
"additive scatter." 

For an explicit scheme, solutions are obtained directly, so that opera- 
tion (2) vectorizes naturally. Operation (3), the application of boundary 
conditions, is an intrinsically scalar operation and difficult to vectorize. 
However, use of bit vectors to flag boundary nodes and use of VPS 32 FORTRAN 
supplied routines enables full vectorization of this operation. 

A program flow chart for the Taylor-Galerkin algorithm is shown in Fig. 
4. 


CONCLUDING REMARKS 

A Taylor-Galerkin finite element solution algorithm for transient non- 
linear thermal-structural analyses of large, complex structural problems 
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subjected to rapidly applied thermal-structural loads is described. The 
two-step Taylor-Galerkin algorithm is an application of an algorithm re- 
cently developed for problems in compressible fluid dynamics. The element 
integrals that appear in the algorithm can be evaluated in closed form for 
two and three dimensional elements. Numerical calculations show that compu- 
tational times are reduced significantly by the closed form integral evalua- 
tion. The algorithm has been implemented on the NASA Langley VPS 32 vector 
computer with special programming strategies to yield very high computa- 
tional speeds for the solution of large problems. The many desirable attri- 
butes of the algorithm indicate that it will be quite effective for the 
solution of large, nonlinear, transient thermal-structural problems. 
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